clear; clc; close all;

k = 3;
c1 = 1e-9; c2 = 0; c3 = 0;
eta = 0;

h = 0.01;
T = 600;
transient = 550;
N = T / h;
N_trans = transient / h;

mu_list = linspace(-3*pi, 3*pi, 300);
LEs_all = zeros(length(mu_list), 3);

parfor i = 1:length(mu_list)
    mu = mu_list(i);
    x0 = [0; 0; 0];
    LEs = computeLEs(k, c1, c2, c3, mu, eta, h, N, N_trans, x0);
    LEs_all(i, :) = LEs(1:3);
end

figure; hold on; box on;
plot(mu_list, LEs_all(:,1), 'r-', 'LineWidth', 1.5);
plot(mu_list, LEs_all(:,2), 'g-', 'LineWidth', 1.5);
yline(0, 'k--');

xlabel('\mu');
ylabel('Lyapunov Exponents');
xticks(-3*pi : pi : 3*pi);
xticklabels({'-3\pi', '-2\pi', '-\pi', '0', '\pi', '2\pi', '3\pi'});
legend('LE_1 (max)', 'LE_2', 'Location', 'best');
grid on;

%% --- 以下函数同之前 ---
function LEs = computeLEs(k, c1, c2, c3, mu, eta, h, N, N_trans, x0)
    n = length(x0);
    X = x0;
    Q = eye(n);
    sum_LE = zeros(n,1);
    
    for step = 1:N
        X = RK4_step(@(x) reduced_system_ode(x,k,c1,c2,c3,mu,eta), X, h);
        J = jacobian_matrix(X,k,mu,eta);
        Q = Q + h*J*Q;
        [Q,R] = qr(Q);
        
        if step > N_trans
            sum_LE = sum_LE + log(abs(diag(R)));
        end
    end
    
    LEs = sum_LE / (h*(N - N_trans));
    LEs = sort(LEs, 'descend')';
end

function dx = reduced_system_ode(x,k,c1,c2,c3,mu,eta)
    X = x(1); Y = x(2); Z = x(3);
    dx = zeros(3,1);
    
    dx(1) = Y + Z - k * sin(Y + eta) + k * sin(eta) + c1;
    dx(2) = -X + Z + c2;
    dx(3) = -X - Z + k * sin(X + mu) - k * sin(mu) + c3;
end

function J = jacobian_matrix(x,k,mu,eta)
    X = x(1); Y = x(2); Z = x(3);
    J = zeros(3,3);
    
    J(1,1) = 0;
    J(1,2) = 1 - k * cos(Y + eta);
    J(1,3) = 1;
    
    J(2,1) = -1;
    J(2,2) = 0;
    J(2,3) = 1;
    
    J(3,1) = -1 + k * cos(X + mu);
    J(3,2) = 0;
    J(3,3) = -1;
end

function x_next = RK4_step(f,x,h)
    k1 = f(x);
    k2 = f(x + h/2*k1);
    k3 = f(x + h/2*k2);
    k4 = f(x + h*k3);
    x_next = x + h/6*(k1 + 2*k2 + 2*k3 + k4);
end
